rm(list=ls())
graphics.off()

library(readxl)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(nlme)
library(matrixcalc)
library(Biodem)
library(stats)
library(multcomp)
library(hrbrthemes)
library(viridis)
library(gapminder)
library(devtools)
library(stringr)
library(patchwork)
library(gridExtra)
library(cowplot)
library(extrafont)
library(gtable)
library(ggtext)
library(effsize)
library(data.table)
library(psych)
library(psychotools)
library(psychTools)
library(BlandAltmanLeh)
library(base)


######################### START ###########################################

setwd("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3")

data <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SC_mean.xlsx")


#Tables
SC_Hand <- bind_cols(data[,1],data[,3:8])
SC_Neck <- bind_cols(data[,1],data[,10:15])

SC_Hand = na.omit(SC_Hand)
SC_Neck = na.omit(SC_Neck)


colnames(SC_Hand) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")
colnames(SC_Neck) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")

###SC_Hand 
SC_Hand_long <- pivot_longer(SC_Hand, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
SC_Hand_long$div <- rep(div,nrow(SC_Hand)) 

SC_Hand_long = split(SC_Hand_long, SC_Hand_long$div)

SC_Hand_coglong = SC_Hand_long$Mental[,1:3]

SC_Hand_motlong =  SC_Hand_long$Motor[,1:3]

###

plot <- ggplot(SC_Hand_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'SC_Hand',
       y = "SC_Hand") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.SC_Hand.CogEasy = shapiro.test(SC_Hand$aclow)
Norm.SC_Hand.CogMiddle = shapiro.test(SC_Hand$bcchall)
Norm.SC_Hand.CogDifficult = shapiro.test(SC_Hand$cchigh)

Norm.SC_Hand.MotEasy = shapiro.test(SC_Hand$dmlow)
Norm.SC_Hand.MotMiddle = shapiro.test(SC_Hand$emchall)
Norm.SC_Hand.MotDifficult = shapiro.test(SC_Hand$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.SC_Hand.CogEasy$p.value > 0.05 & Norm.SC_Hand.CogMiddle$p.value > 0.05 & Norm.SC_Hand.CogDifficult$p.value > 0.05){
  print(anova_test(data=SC_Hand_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(SC_Hand_coglong$value, SC_Hand_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SC_Hand_coglong, value ~ name|ID))
  pairwise.wilcox.test(SC_Hand_coglong$value, SC_Hand_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.SC_Hand.MotEasy$p.value > 0.05 & Norm.SC_Hand.MotMiddle$p.value > 0.05 & Norm.SC_Hand.MotDifficult$p.value > 0.05){
  print(anova_test(data=SC_Hand_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(SC_Hand_motlong$value, SC_Hand_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SC_Hand_motlong, value ~ name|ID))
  pairwise.wilcox.test(SC_Hand_motlong$value, SC_Hand_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_SC_Hand_Test <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SC_Day1.xlsx")
data_SC_Hand_ReTest <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SC_Day2.xlsx")


## SC_Hand

## Cognitive

### Easy
data_SC_Hand_Cog_easy = data.table(data_SC_Hand_Test$SC_Cog_Hand_Easy, data_SC_Hand_ReTest$SC_Cog_Hand_Easy)
#data_SC_Hand_easy_RRn = data_SC_Hand_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Cog_easy =ICC(data_SC_Hand_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Cog_easy<- bland.altman.plot(data_SC_Hand_Test$SC_Cog_Hand_Easy,data_SC_Hand_ReTest$SC_Cog_Hand_Easy, main="Test-Retest SC_Hand Easy", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Cog_easy)

###Middle
data_SC_Hand_Cog_middle = data.table(data_SC_Hand_Test$SC_Cog_Hand_Optimal, data_SC_Hand_ReTest$SC_Cog_Hand_Optimal)
#data_SC_Hand_Cog_middle = data_SC_Hand_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Cog_middle =ICC(data_SC_Hand_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Cog_middle <-bland.altman.plot(data_SC_Hand_Test$SC_Cog_Hand_Optimal,data_SC_Hand_ReTest$SC_Cog_Hand_Optimal, main="Test-Retest SC_Hand Optimal", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Cog_middle)


###Difficult
data_SC_Hand_Cog_difficult = data.table(data_SC_Hand_Test$SC_Cog_Hand_Difficult, data_SC_Hand_ReTest$SC_Cog_Hand_Difficult)
#data_SC_Hand_Cog_difficult = data_SC_Hand_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Cog_difficult =ICC(data_SC_Hand_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Cog_difficult<- bland.altman.plot(data_SC_Hand_Test$SC_Cog_Hand_Difficult,data_SC_Hand_ReTest$SC_Cog_Hand_Difficult, main="Test-Retest SC_Hand Difficult", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Cog_difficult)

## Motor

### Easy
data_SC_Hand_Mot_easy = data.table(data_SC_Hand_Test$SC_Mot_Hand_Easy, data_SC_Hand_ReTest$SC_Mot_Hand_Easy)
#data_SC_Hand_easy_RRn = data_SC_Hand_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Mot_easy =ICC(data_SC_Hand_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Mot_easy<- bland.altman.plot(data_SC_Hand_Test$SC_Mot_Hand_Easy,data_SC_Hand_ReTest$SC_Mot_Hand_Easy, main="Test-Retest SC_Hand Easy", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Mot_easy)

###Middle
data_SC_Hand_Mot_middle = data.table(data_SC_Hand_Test$SC_Mot_Hand_Optimal, data_SC_Hand_ReTest$SC_Mot_Hand_Optimal)
#data_SC_Hand_Mot_middle = data_SC_Hand_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Mot_middle =ICC(data_SC_Hand_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Mot_middle <-bland.altman.plot(data_SC_Hand_Test$SC_Mot_Hand_Optimal,data_SC_Hand_ReTest$SC_Mot_Hand_Optimal, main="Test-Retest SC_Hand Optimal", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Mot_middle)


###Difficult
data_SC_Hand_Mot_difficult = data.table(data_SC_Hand_Test$SC_Mot_Hand_Difficult, data_SC_Hand_ReTest$SC_Mot_Hand_Difficult)
#data_SC_Hand_Mot_difficult = data_SC_Hand_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Hand_Mot_difficult =ICC(data_SC_Hand_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Hand_Mot_difficult<- bland.altman.plot(data_SC_Hand_Test$SC_Mot_Hand_Difficult,data_SC_Hand_ReTest$SC_Mot_Hand_Difficult, main="Test-Retest SC_Hand Difficult", xlab="Means", ylab="differences")

print(ICC_SC_Hand_Mot_difficult)


###################################################################################################################################################


###SC_Neck 
SC_Neck_long <- pivot_longer(SC_Neck, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
SC_Neck_long$div <- rep(div,nrow(SC_Neck)) 

SC_Neck_long = split(SC_Neck_long, SC_Neck_long$div)

SC_Neck_coglong = SC_Neck_long$Mental[,1:3]

SC_Neck_motlong =  SC_Neck_long$Motor[,1:3]

###

plot <- ggplot(SC_Neck_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'SC_Neck',
       y = "SC_Neck") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.SC_Neck.CogEasy = shapiro.test(SC_Neck$aclow)
Norm.SC_Neck.CogMiddle = shapiro.test(SC_Neck$bcchall)
Norm.SC_Neck.CogDifficult = shapiro.test(SC_Neck$cchigh)

Norm.SC_Neck.MotEasy = shapiro.test(SC_Neck$dmlow)
Norm.SC_Neck.MotMiddle = shapiro.test(SC_Neck$emchall)
Norm.SC_Neck.MotDifficult = shapiro.test(SC_Neck$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.SC_Neck.CogEasy$p.value > 0.05 & Norm.SC_Neck.CogMiddle$p.value > 0.05 & Norm.SC_Neck.CogDifficult$p.value > 0.05){
  print(anova_test(data=SC_Neck_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(SC_Neck_coglong$value, SC_Neck_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SC_Neck_coglong, value ~ name|ID))
  pairwise.wilcox.test(SC_Neck_coglong$value, SC_Neck_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.SC_Neck.MotEasy$p.value > 0.05 & Norm.SC_Neck.MotMiddle$p.value > 0.05 & Norm.SC_Neck.MotDifficult$p.value > 0.05){
  print(anova_test(data=SC_Neck_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(SC_Neck_motlong$value, SC_Neck_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SC_Neck_motlong, value ~ name|ID))
  pairwise.wilcox.test(SC_Neck_motlong$value, SC_Neck_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_SC_Neck_Test <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SC_Day1.xlsx")
data_SC_Neck_ReTest <- read_xlsx("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SC_Day2.xlsx")


## SC_Neck

## Cognitive

### Easy
data_SC_Neck_Cog_easy = data.table(data_SC_Neck_Test$SC_Cog_Neck_Easy, data_SC_Neck_ReTest$SC_Cog_Neck_Easy)
#data_SC_Neck_easy_RRn = data_SC_Neck_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Cog_easy =ICC(data_SC_Neck_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Cog_easy<- bland.altman.plot(data_SC_Neck_Test$SC_Cog_Neck_Easy,data_SC_Neck_ReTest$SC_Cog_Neck_Easy, main="Test-Retest SC_Neck Easy", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Cog_easy)

###Middle
data_SC_Neck_Cog_middle = data.table(data_SC_Neck_Test$SC_Cog_Neck_Optimal, data_SC_Neck_ReTest$SC_Cog_Neck_Optimal)
#data_SC_Neck_Cog_middle = data_SC_Neck_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Cog_middle =ICC(data_SC_Neck_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Cog_middle <-bland.altman.plot(data_SC_Neck_Test$SC_Cog_Neck_Optimal,data_SC_Neck_ReTest$SC_Cog_Neck_Optimal, main="Test-Retest SC_Neck Optimal", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Cog_middle)


###Difficult
data_SC_Neck_Cog_difficult = data.table(data_SC_Neck_Test$SC_Cog_Neck_Difficult, data_SC_Neck_ReTest$SC_Cog_Neck_Difficult)
#data_SC_Neck_Cog_difficult = data_SC_Neck_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Cog_difficult =ICC(data_SC_Neck_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Cog_difficult<- bland.altman.plot(data_SC_Neck_Test$SC_Cog_Neck_Difficult,data_SC_Neck_ReTest$SC_Cog_Neck_Difficult, main="Test-Retest SC_Neck Difficult", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Cog_difficult)

## Motor

### Easy
data_SC_Neck_Mot_easy = data.table(data_SC_Neck_Test$SC_Mot_Neck_Easy, data_SC_Neck_ReTest$SC_Mot_Neck_Easy)
#data_SC_Neck_easy_RRn = data_SC_Neck_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Mot_easy =ICC(data_SC_Neck_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Mot_easy<- bland.altman.plot(data_SC_Neck_Test$SC_Mot_Neck_Easy,data_SC_Neck_ReTest$SC_Mot_Neck_Easy, main="Test-Retest SC_Neck Easy", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Mot_easy)

###Middle
data_SC_Neck_Mot_middle = data.table(data_SC_Neck_Test$SC_Mot_Neck_Optimal, data_SC_Neck_ReTest$SC_Mot_Neck_Optimal)
#data_SC_Neck_Mot_middle = data_SC_Neck_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Mot_middle =ICC(data_SC_Neck_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Mot_middle <-bland.altman.plot(data_SC_Neck_Test$SC_Mot_Neck_Optimal,data_SC_Neck_ReTest$SC_Mot_Neck_Optimal, main="Test-Retest SC_Neck Optimal", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Mot_middle)


###Difficult
data_SC_Neck_Mot_difficult = data.table(data_SC_Neck_Test$SC_Mot_Neck_Difficult, data_SC_Neck_ReTest$SC_Mot_Neck_Difficult)
#data_SC_Neck_Mot_difficult = data_SC_Neck_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SC_Neck_Mot_difficult =ICC(data_SC_Neck_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SC_Neck_Mot_difficult<- bland.altman.plot(data_SC_Neck_Test$SC_Mot_Neck_Difficult,data_SC_Neck_ReTest$SC_Mot_Neck_Difficult, main="Test-Retest SC_Neck Difficult", xlab="Means", ylab="differences")

print(ICC_SC_Neck_Mot_difficult)

